Metode Spasiotemporal Cokriging untuk Penaksiran Konsentrasi \(PM_{2.5}\) di Provinsi Daerah Khusus Ibukota Jakarta

Pendahuluan

Polusi udara merupakan masalah lingkungan yang signifikan di kota-kota padat penduduk, terutama yang menjadi pusat kegiatan ekonomi dan industri, seperti Jakarta, ibu kota Indonesia. Polusi udara di Jakarta sebagian besar berasal dari konsumsi energi, penggunaan bahan bakar fosil, aktivitas rumah tangga, dan proses industri. Di antara polutan udara, partikel materi halus 2.5 (\(PM_{2.5}\)) menjadi ancaman serius bagi kesehatan manusia karena ukurannya yang kecil (≤ 2,5 mikrometer, disingkat µm), yang memungkinkan untuk masuk jauh ke dalam sistem pernapasan. Menyadari bahayanya, Pemerintah Indonesia telah menetapkan ambang batas sehat untuk \(PM_{2.5}\) sebesar 50 µg/m³ (mikrogram per meter kubik). Informasi yang dapat diandalkan tentang konsentrasi \(PM_{2.5}\) sangat penting bagi pembuat kebijakan untuk menyusun strategi mitigasi dan bagi masyarakat untuk melindungi diri dari paparan.

Meskipun penting, pemantauan konsentrasi \(PM_{2.5}\) sering kali terbatas oleh ketersediaan peralatan pemantauan, yang mengakibatkan kekosongan data di berbagai wilayah. Untuk mengatasi hal ini, teknik interpolasi digunakan untuk memperkirakan konsentrasi polutan di lokasi yang tidak terpantau. Metode interpolasi yang umum digunakan meliputi Inverse Distance Weighting (IDW) dan kriging.

IDW memperkirakan nilai di titik yang tidak diketahui hanya berdasarkan jaraknya dari lokasi yang diketahui, dengan memberikan bobot lebih tinggi pada titik yang lebih dekat. Namun, IDW tidak memperhitungkan ketergantungan spasial atau pola variabilitas yang mendasarinya, dan tidak dapat memasukkan variabel prediktor tambahan, yang membatasi kemampuannya untuk memberikan estimasi akurat dalam sistem lingkungan yang kompleks.

Kriging meningkatkan IDW dengan memanfaatkan ketergantungan spasial melalui pemodelan semivariogram, yang menangkap hubungan antara variabilitas dan jarak di antara lokasi yang diamati. Hal ini menghasilkan prediksi yang lebih akurat dibandingkan IDW, terutama di wilayah dengan autokorelasi spasial yang signifikan. Namun, seperti IDW, kriging tidak memungkinkan dimasukkannya variabel prediktor dan tidak secara eksplisit mempertimbangkan ketergantungan temporal.

Namun, polusi udara dipengaruhi oleh kombinasi faktor spasial dan temporal, serta variabel sekunder seperti curah hujan, kelembapan, kecepatan angin, dan konsentrasi \(NO_2\). Untuk mengatasi keterbatasan ini, cokriging spasiotemporal digunakan karena mengintegrasikan ketergantungan spasial dan temporal ke dalam model serta memungkinkan dimasukkannya variabel sekunder sebagai prediktor. Pendekatan ini memungkinkan analisis yang lebih komprehensif terhadap pola distribusi polutan dan meningkatkan akurasi prediksi dengan memasukkan pengaruh faktor lingkungan tambahan bersama konsentrasi \(PM_{2.5}\).

Dalam konteks ini, mengintegrasikan ketergantungan spasial-temporal dan variabel lingkungan tambahan ke dalam model cokriging diharapkan dapat meningkatkan akurasi estimasi konsentrasi \(PM_{2.5}\). Penelitian ini bertujuan untuk mengembangkan model cokriging spasiotemporal untuk memberikan prediksi yang lebih akurat tentang konsentrasi \(PM_{2.5}\) di Jakarta. Dengan memasukkan curah hujan, kelembapan, kecepatan angin, dan konsentrasi \(NO_2\) sebagai variabel sekunder, model yang diusulkan bertujuan untuk menangkap variasi spasial dan temporal secara lebih efektif, sehingga menghasilkan estimasi yang lebih akurat yang mencerminkan dinamika polusi.

Load package

library(gstat); library(sf); library(spacetime); library(sp);library(ggplot2)

## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE

library(plotly); library(xts); library(RColorBrewer); library(openxlsx)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(raster); library(geojson); library(patchwork); library(ggplotify)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:plotly':
## 
##     select

## 
## Attaching package: 'geojson'

## The following object is masked from 'package:graphics':
## 
##     polygon

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:raster':
## 
##     area

library(moments); library(tidyr); library(dplyr); library(psych)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:raster':
## 
##     extract

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:xts':
## 
##     first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## 
## Attaching package: 'psych'

## The following object is masked from 'package:raster':
## 
##     distance

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

Statistik Deskriptif

Pertama, analisis statistik deskriptif dilakukan pada variabel-variabel penelitian.

# INPUT DATA
# PM2.5
pm25 = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx", 
                 sheet = "PM25")
pm25$time = ISOdate(pm25$Tahun, pm25$Bulan, 1, tz = "GMT")

stations = 3:31
pm25_mat = as.matrix(pm25[stations])
# Curah Hujan
presipitasi = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx", 
                        sheet = "Presipitasi")
presipitasi$time = ISOdate(presipitasi$Tahun, presipitasi$Bulan, 1, tz = "GMT")
stations = 3:31
presipitasi_mat = as.matrix(presipitasi[stations])
# Humidity
kelembaban = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx", 
                       sheet = "Kelembaban")  
kelembaban$time = ISOdate(kelembaban$Tahun, kelembaban$Bulan, 1, tz = "GMT")
stations = 3:31
kelembaban_mat = as.matrix(kelembaban[stations])
# Data Kecepatan Angin
angin = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx", 
                  sheet = "Kecepatan Angin")
angin$time = ISOdate(angin$Tahun, angin$Bulan, 1, tz = "GMT")
stations = 3:31
angin_mat = as.matrix(angin[stations])
# Data NO2
no2 = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx", 
                sheet = "NO2")
no2$time = ISOdate(no2$Tahun, no2$Bulan, 1, tz = "GMT")
stations = 3:31
no2_mat = as.matrix(no2[stations])

# MEMBUAT STFDF CLASS OBJECT
tesis.loc = read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx",
                      sheet = "Koordinat")
coordinates(tesis.loc) = ~Longitude+Latitude
proj4string(tesis.loc) = "+proj=longlat +datum=WGS84"
pts = coordinates(tesis.loc)
rownames(pts) = tesis.loc$Sensor
pts = SpatialPoints(pts, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84"))
utm48s = "+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"
pts.sfc = st_transform(st_as_sfc(pts), utm48s)
pts = as(pts.sfc, "Spatial")
# membuat spatial object dari data masing-masing variabel
pm25.data = stConstruct(pm25_mat, space = list(values = 1:ncol(pm25_mat)),
                        time = pm25$time, SpatialObj = pts, interval = TRUE)
presipitasi.data = stConstruct(presipitasi_mat, space = list(values = 1:ncol(presipitasi_mat)),
                               time = presipitasi$time, SpatialObj = pts, interval = TRUE)
kelembaban.data = stConstruct(kelembaban_mat, space = list(values = 1:ncol(kelembaban_mat)),
                              time = kelembaban$time, SpatialObj = pts, interval = TRUE)
angin.data = stConstruct(angin_mat, space = list(values = 1:ncol(angin_mat)),
                         time = angin$time, SpatialObj = pts, interval = TRUE)
no2.data = stConstruct(no2_mat, space = list(values = 1:ncol(no2_mat)),
                       time = no2$time, SpatialObj = pts, interval = TRUE)
# STFDF Object untuk gabungan
combined_data <- cbind(as.data.frame(pm25.data)[,-7],
                       PM25 = as.data.frame(pm25.data)[,7],
                       Presipitasi = as.data.frame(presipitasi.data)[,7],
                       Kelembaban = as.data.frame(kelembaban.data)[,7],
                       Angin = as.data.frame(angin.data)[,7],
                       NO2 = as.data.frame(no2.data)[,7])

descriptive_stats <- function(x) {
  c(
    Min = min(x, na.rm = TRUE),
    Max = max(x, na.rm = TRUE),
    Mean = mean(x, na.rm = TRUE),
    `Std Dev` = sd(x, na.rm = TRUE)
  )
}

# Menghitung statistik deskriptif dan membuat hasil dalam format panjang.
results_long <- combined_data[,6:11] %>%
  pivot_longer(cols = -timeIndex, names_to = "Variable", values_to = "Value") %>%
  group_by(Variable) %>%
  summarise(
    Minimum = min(Value, na.rm = TRUE),
    Maximum = max(Value, na.rm = TRUE),
    Mean = mean(Value, na.rm = TRUE),
    `Std Dev` = sd(Value, na.rm = TRUE),
    .groups = "drop"
  )

# Menampilkan hasil
print(results_long)

## # A tibble: 5 × 5
##   Variable    Minimum Maximum   Mean `Std Dev`
##   <chr>         <dbl>   <dbl>  <dbl>     <dbl>
## 1 Angin          2.14    6.22   3.79     0.870
## 2 Kelembaban    67.1    89.8   78.4      4.83 
## 3 NO2            3.9    30.0   19.4      5.83 
## 4 PM25          14.6    68.8   39.6     11.7  
## 5 Presipitasi    0     559.   140.     138.

Konsentrasi \(PM_{2.5}\) menunjukkan variasi sepanjang tahun, yang ditunjukkan oleh standar deviasi sebesar 11,69 µg/m³. Nilai ini menunjukkan bahwa konsentrasi \(PM_{2.5}\) menyimpang sebesar 11,69 µg/m³ dari rata-rata. Konsentrasi \(PM_{2.5}\) terendah terjadi pada bulan Februari, dengan nilai 14,56 µg/m³, sedangkan konsentrasi maksimum tercatat pada bulan Agustus, mencapai 68,83 µg/m³. Konsentrasi \(NO_{2}\) terendah adalah 3,9 µg/m³, yang terjadi pada bulan Desember, sedangkan konsentrasi tertinggi, juga tercatat pada bulan Desember tetapi di lokasi yang berbeda, adalah 29,97 µg/m³. Curah hujan menunjukkan variabilitas yang lebih tinggi dibandingkan dengan variabel lainnya, sebagaimana tercermin dalam standar deviasi sebesar 137,78 mm. Curah hujan minimum di Jakarta selama periode 2023 adalah 0 mm, yang menunjukkan tidak ada hujan sama sekali dalam satu bulan penuh, sedangkan curah hujan maksimum mencapai 558,83 mm. Variabel kelembapan udara dan kecepatan angin menunjukkan nilai yang relatif homogen di seluruh lokasi dan periode waktu, dengan standar deviasi masing-masing sebesar 4,82% dan 0,87 m/s.

Analisis Korelasi

Terdapat korelasi negatif antara curah hujan, kelembapan udara, dan kecepatan angin terhadap konsentrasi \(PM_{2.5}\). Ketika kelembapan udara, kecepatan angin, dan curah hujan meningkat, konsentrasi \(PM_{2.5}\) cenderung menurun. Sebaliknya, \(NO_{2}\) memiliki korelasi positif dengan \(PM_{2.5}\). Ini berarti bahwa konsentrasi \(PM_{2.5}\) juga akan meningkat seiring dengan meningkatnya \(NO_{2}\).

corr.test(combined_data[,c(7,11,8:10)], method = 'pearson', adjust = 'none')

## Call:corr.test(x = combined_data[, c(7, 11, 8:10)], method = "pearson", 
##     adjust = "none")
## Correlation matrix 
##              PM25   NO2 Presipitasi Kelembaban Angin
## PM25         1.00  0.31       -0.58      -0.13 -0.75
## NO2          0.31  1.00       -0.34      -0.43 -0.19
## Presipitasi -0.58 -0.34        1.00       0.73  0.54
## Kelembaban  -0.13 -0.43        0.73       1.00  0.04
## Angin       -0.75 -0.19        0.54       0.04  1.00
## Sample Size 
## [1] 348
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##             PM25 NO2 Presipitasi Kelembaban Angin
## PM25        0.00   0           0       0.02  0.00
## NO2         0.00   0           0       0.00  0.00
## Presipitasi 0.00   0           0       0.00  0.00
## Kelembaban  0.02   0           0       0.00  0.47
## Angin       0.00   0           0       0.47  0.00
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Data Splitting

Dari total sampel sebanyak 29 titik, data akan dibagi menjadi 80% data pelatihan (23 titik) dan 20% data pengujian (6 titik).

# Baca data lokasi stasiun dari file Excel
sensor_loc <- read.xlsx("D:/TESIS/DATA/Data Fix/DataRequest_UNPAD_FarizaAP_Jakarta2023.xlsx", 
                        sheet = 'Koordinat')
# RANDOM SAMPLE
sample_index = unique(combined_data$sp.ID)
#set.seed(123) # untuk konsistensi
#test_indices <- sample(sample_index, size = round(0.2 * length(sample_index),0))
#saveRDS(test_indices, 
#        "D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
test_samples = readRDS("D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
train_samples = as.numeric(setdiff(sample_index, test_samples))

# Convert the data frame to a spatial points data frame
stasiun_test = st_as_sf(sensor_loc[sensor_loc$Sensor %in% c("Cibubur", "Cipinang Besar",
                                                      "Duri Utara", "Marunda",
                                                      "Rempoa Permai", "Semanan"),], 
                         coords = c("Longitude", "Latitude"), crs = 4326)
stasiun_train = st_as_sf(sensor_loc[!(sensor_loc$Sensor %in% c("Cibubur", "Cipinang Besar",
                                                            "Duri Utara", "Marunda",
                                                            "Rempoa Permai", "Semanan")),],
                        coords = c("Longitude", "Latitude"), crs = 4326)

# Load the Jakarta boundary shapefile (adjust the path to your shapefile)
province_boundary = st_read("D:/TESIS/DATA/SHP/gadm41_IDN_2.shp")

## Reading layer `gadm41_IDN_2' from data source `D:\TESIS\DATA\SHP\gadm41_IDN_2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS:  WGS 84

# Inspect the structure of the province_boundary
str(province_boundary)

## Classes 'sf' and 'data.frame':   502 obs. of  14 variables:
##  $ GID_2    : chr  "IDN.1.2_1" "IDN.1.1_1" "IDN.1.3_1" "IDN.1.4_1" ...
##  $ GID_0    : chr  "IDN" "IDN" "IDN" "IDN" ...
##  $ COUNTRY  : chr  "Indonesia" "Indonesia" "Indonesia" "Indonesia" ...
##  $ GID_1    : chr  "IDN.1_1" "IDN.1_1" "IDN.1_1" "IDN.1_1" ...
##  $ NAME_1   : chr  "Aceh" "Aceh" "Aceh" "Aceh" ...
##  $ NL_NAME_1: chr  "NA" "NA" "NA" "NA" ...
##  $ NAME_2   : chr  "Aceh Barat" "Aceh Barat Daya" "Aceh Besar" "Aceh Jaya" ...
##  $ VARNAME_2: chr  "NA" "NA" "NA" "NA" ...
##  $ NL_NAME_2: chr  "NA" "NA" "NA" "NA" ...
##  $ TYPE_2   : chr  "Kabupaten" "Kabupaten" "Kabupaten" "Kabupaten" ...
##  $ ENGTYPE_2: chr  "Regency" "Regency" "Regency" "Regency" ...
##  $ CC_2     : chr  "1107" "1112" "1108" "1116" ...
##  $ HASC_2   : chr  "ID.AC.AB" "ID.AC.AD" "ID.AC.AR" "ID.AC.AJ" ...
##  $ geometry :sfc_MULTIPOLYGON of length 502; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:929, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:13] "GID_2" "GID_0" "COUNTRY" "GID_1" ...

jkt_boundary <- province_boundary %>% 
  filter(NAME_1 == "Jakarta Raya")
#jkt_boundary[-6,"NAME_EN"] = c('West Jakarta', 'Central Jakarta', 'South Jakarta',
#                               'East Jakarta', 'North Jakarta')
jkt_boundary[-6,"NAME_EN"] = c('Jakarta Barat', 'Jakarta Pusat', 'Jakarta Selatan',
                               'Jakarta Timur', 'Jakarta Utara')
# Extract the boundary lines for visualization
boundary_lines <- st_cast(jkt_boundary, "MULTILINESTRING")

# Plot the points and the filtered Jakarta boundary with legend title as "SPKU"
ggplot() +
  geom_sf(data = jkt_boundary[-6,], fill = "lightyellow", color = "brown") +
  geom_sf_text(data = jkt_boundary[-6,], aes(label = NAME_EN), size = 4, 
               nudge_y = c(0.01, 0.01, 0.02, 0.02, 0.01), 
               nudge_x = c(0, 0.02, 0, 0.02, 0), color = "black") +  
  geom_sf(data = stasiun_train, aes(color = "Latih"), size = 2) +  # Tambahkan label legenda
  geom_sf(data = stasiun_test, aes(color = "Uji"), size = 2) +    # Tambahkan label legenda
  scale_color_manual(name = "Data", values = c("Latih" = "blue", "Uji" = "orangered")) + # Skala warna
  theme_minimal() +
  labs(x = "Longitude",
       y = "Latitude") +
  theme(panel.grid.major = element_blank(),  # Ubah grid besar jadi putih
        panel.grid.minor = element_blank(),  # Ubah grid kecil jadi putih
        legend.title = element_text(size = 12),  # Ukuran judul legenda
        legend.text = element_text(size = 10))   # Ukuran teks legenda

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Gambar tersebut menunjukkan pembagian data yang digunakan sebagai data latih dan data uji. Observasi yang ditunjukkan oleh titik berwarna biru menunjukkan data latih, sedangkan titik berwarna merah menunjukkan data uji. Berikut ini adalah proses membagi data dan mengubah data ke dalam format STFDF

# Konversi kolom waktu ke format POSIXct
Sys.setenv(TZ = "GMT")
time = as.POSIXct(unique(combined_data$time), tz = "GMT")

# RANDOM SAMPLE
sample_index = unique(combined_data$sp.ID)
#set.seed(123) # untuk konsistensi
#test_indices <- sample(sample_index, size = round(0.2 * length(sample_index),0))
#saveRDS(test_indices, 
#        "D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")
test_samples = readRDS("D:/TESIS/DATA/Data Fix/indeks wilayah test.RDS")

# Membuat objek STFDF DATA TRAIN
train_samples = as.numeric(setdiff(sample_index, test_samples))
train_data = combined_data[combined_data$sp.ID %in% train_samples, ]
coords_train = unique(train_data[, c("coords.x1","coords.x2")])
pts_train = coordinates(coords_train)
pts_train = SpatialPoints(pts_train,
                          CRS("+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"))
utm48s = "+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"
pts_train.sfc = st_transform(st_as_sfc(pts_train), utm48s)
pts_train = as(pts_train.sfc, "Spatial")
stfdf.train = STFDF(sp = pts_train, time = time, data = train_data[,7:11])
stfdf.train

## An object of class "STFDF"
## Slot "data":
##         PM25 Presipitasi Kelembaban Angin      NO2
## 1   25.26000      169.67   80.50000  4.87  9.37000
## 2   29.50000      185.37   82.27000  4.64  7.97000
## 3   24.23000      165.96   80.07000  4.92  9.71000
## 5   31.47000      192.90   83.10000  4.53  7.31000
## 6   23.65000      163.89   79.83000  4.95  9.90000
## 8   28.92000      183.18   82.03000  4.67  8.16000
## 10  21.32000      155.70   78.85000  5.08 10.67000
## 11  21.18000      155.22   78.79000  5.09 10.71000
## 12  23.93667      164.91   79.95000  4.94  9.80000
## 13  26.77000      175.18   81.13000  4.78  8.87000
## 14  20.47000      152.77   78.50000  5.13 21.60000
## 15  21.39000      155.95   78.88000  5.08 10.64000
## 16  20.42000      152.60   78.48000  5.13 10.96000
## 17  28.30000      180.86   81.77000  4.70  8.36000
## 19  22.67000      160.42   79.42000  5.01 10.22000
## 20  17.01000      141.11   77.05000  5.31 12.09000
## 21  29.31000      184.65   82.19000  4.65  8.03000
## 22  22.77000      160.77   79.46000  5.00 10.19000
## 23  22.07000      158.32   79.17000  5.04 10.42000
## 26  18.78000      147.02   77.79000  5.22 11.51000
## 27  21.62000      169.57   85.69774  5.40 10.57000
## 28  22.68000      160.46   79.42000  5.01 10.22000
## 29  22.80000      160.88   79.47000  5.00 10.18000
## 30  20.87000      502.43   84.09000  5.88 14.36000
## 31  23.76000      520.60   85.30000  5.72 13.41000
## 32  21.88000      508.74   84.52000  5.82 14.03000
## 34  23.49000      518.89   85.19000  5.73 13.50000
## 35  20.35000      499.20   83.88000  5.90 14.54000
## 37  24.85000      527.54   85.76000  5.66 13.05000
## 39  19.31333      492.78   83.44000  5.96 14.88000
## 40  21.85000      508.56   84.50000  5.82 14.04000
## 41  21.23333      504.70   84.25000  5.86 14.24000
## 42  23.01000      515.85   84.99000  5.76 13.66000
## 43  17.78000      483.36   82.80000  6.04 18.64286
## 44  20.34000      499.13   83.87000  5.90 14.54000
## 45  17.92000      484.22   82.86000  6.04 15.34000
## 46  23.57000      519.40   85.22000  5.73 13.47000
## 48  20.51000      500.19   83.94000  5.89 14.48000
## 49  14.56000      463.88   81.46000  6.22 16.45000
## 50  26.45000      537.80   86.43000  5.57 12.52000
## 51  18.59000      488.32   83.14000  6.00 15.12000
## 52  19.28000      492.57   83.43000  5.96 14.89000
## 55  20.15000      497.95   83.79000  5.91 14.60000
## 56  18.82000      442.02   86.65000  5.60 15.04000
## 57  20.23000      498.45   83.83000  5.91 14.58000
## 58  21.18000      504.36   84.22000  5.86 14.26000
## 59  39.89000      285.12   86.16000  3.54 20.24000
## 60  38.60000      279.10   85.62000  3.61 20.66000
## 61  31.84000      248.58   82.79000  3.97 22.89000
## 63  42.27000      296.40   87.16000  3.41 19.45000
## 64  32.73000      252.50   83.17000  3.92 22.60000
## 66  40.03000      285.78   86.22000  3.53 20.19000
## 68  30.09667      241.00   82.07000  4.07 23.47000
## 69  32.92000      253.34   83.25000  3.91 22.54000
## 70  29.21000      237.19   81.69000  4.12 23.76000
## 71  39.38000      282.73   85.95000  3.56 20.40000
## 72  31.91000      248.89   82.82000  3.97 24.80645
## 73  31.48000      247.01   82.64000  3.99 23.01000
## 74  28.74000      235.18   81.50000  4.14 23.92000
## 75  37.90000      275.85   85.33000  3.64 20.89000
## 77  29.69000      239.25   81.90000  4.09 23.60000
## 78  23.83000      214.72   79.45000  4.41 25.54000
## 79  37.80000      275.39   85.29000  3.65 20.93000
## 80  31.27000      246.09   82.56000  4.00 23.08000
## 81  26.89000      227.36   80.72000  4.24 24.53000
## 84  33.18000      254.49   83.35000  3.90 22.45000
## 85  27.82000      231.10   84.35871  3.40 24.22000
## 86  28.18000      232.80   81.26000  4.17 24.10000
## 87  28.65000      234.80   81.46000  4.15 23.95000
## 88  32.54000      182.31   80.93000  3.47 15.78000
## 89  34.40000      189.36   81.71000  3.37 15.17000
## 90  25.49000      156.81   77.99000  3.85 18.11000
## 92  37.67000      202.08   83.08000  3.19 14.09000
## 93  28.02000      165.74   79.04000  3.72 17.28000
## 95  34.63000      190.24   81.81000  3.36 15.09000
## 97  24.70000      154.07   77.66000  3.90 18.37000
## 98  26.77000      161.30   78.52000  3.78 17.69000
## 99  25.49000      156.81   77.99000  3.85 18.11000
## 100 33.20000      184.80   81.21000  3.44 15.56000
## 101 28.21000      166.42   79.12000  3.71 18.86667
## 102 23.85000      151.15   77.30000  3.94 18.65000
## 103 25.16000      155.66   77.85000  3.87 18.22000
## 104 32.02000      180.37   80.72000  3.50 15.95000
## 106 25.87000      158.14   78.15000  3.83 17.99000
## 107 22.79000      147.55   76.86000  4.00 19.00000
## 108 33.57000      186.20   81.37000  3.42 15.44000
## 109 28.30000      166.74   79.16000  3.70 17.18000
## 110 23.42000      149.69   77.12000  3.97 18.79000
## 113 29.89000      172.50   79.83000  3.61 16.66000
## 114 26.12000      114.66   83.01833  4.50 17.90000
## 115 25.98000      158.52   78.19000  3.83 17.95000
## 116 26.13000      159.05   78.25000  3.82 17.90000
## 117 55.43000      164.05   84.00000  2.72 20.84000
## 118 58.48000      175.09   85.27000  2.56 19.84000
## 119 36.25000      102.87   75.98000  3.76 27.18000
## 121 60.11000      181.14   85.95000  2.47 19.30000
## 122 44.82000      128.45   79.56000  3.30 24.35000
## 124 53.24000      156.35   83.08000  2.84 21.57000
## 126 49.44000      143.42   81.49000  3.05 22.82000
## 127 46.85000      134.92   80.41000  3.19 23.68000
## 128 39.72000      112.88   77.43000  3.57 26.03000
## 129 53.34000      156.69   83.12000  2.83 21.53000
## 130 36.65000      104.00   76.14000  3.74 22.56667
## 131 42.10000      120.02   78.42000  3.44 25.25000
## 132 41.65000      118.66   78.23000  3.47 25.40000
## 133 50.75000      147.81   82.04000  2.98 22.39000
## 135 40.28000      114.54   77.66000  3.54 25.85000
## 136 37.94000      107.69   76.68000  3.67 26.62000
## 137 49.91000      144.99   81.69000  3.02 22.67000
## 138 46.94000      135.22   80.45000  3.18 23.65000
## 139 34.64000       98.38   75.30000  3.85 27.71000
## 142 46.32000      133.22   80.19000  3.22 23.85000
## 143 42.95000      172.48   81.13839  3.80 24.97000
## 144 32.01000       91.26   74.20000  3.99 28.58000
## 145 40.56000      115.38   77.78000  3.53 25.76000
## 146 53.91000      161.11   84.71000  2.56 17.82000
## 147 49.31000      145.29   82.79000  2.81 19.34000
## 148 38.25000      110.59   78.17000  3.41 22.99000
## 150 53.99000      161.39   84.75000  2.56 17.79000
## 151 41.67000      120.81   79.60000  3.22 21.86000
## 153 51.57000      152.96   83.74000  2.69 18.59000
## 155 42.84000      124.42   80.09000  3.16 21.48000
## 156 42.24000      122.56   79.83000  3.19 21.68000
## 157 40.42000      117.02   79.07000  3.29 22.28000
## 158 49.76000      146.80   82.98000  2.78 19.19000
## 159 43.93000      127.82   80.54000  3.10 22.13333
## 160 41.61000      120.63   79.57000  3.23 21.88000
## 161 40.73000      117.96   79.20000  3.27 22.17000
## 162 48.47000      142.48   82.44000  2.85 19.62000
## 164 39.09000      113.06   78.52000  3.36 22.72000
## 165 35.60000      102.98   77.06000  3.55 23.87000
## 166 50.86000      150.53   83.44000  2.73 18.83000
## 167 45.93000      134.18   81.38000  2.99 20.46000
## 168 34.65000      100.32   76.66000  3.60 24.18000
## 171 42.48000      123.30   79.94000  3.18 21.60000
## 172 43.32000      107.15   82.91300  3.70 21.32000
## 173 36.72000      106.16   77.53000  3.49 23.50000
## 174 40.03000      115.85   78.91000  3.31 22.41000
## 175 56.96000       27.28   80.30000  3.02 20.62000
## 176 51.91000       20.44   78.19000  3.29 22.28000
## 177 38.28000        6.90   72.49000  4.03 26.79000
## 179 67.43000       44.59   84.68000  2.45 17.16000
## 180 47.15000       14.89   76.20000  3.55 23.86000
## 182 55.72000       25.51   79.79000  3.09 21.03000
## 184 41.39667        9.36   73.80000  3.86 25.76000
## 185 48.07000       15.90   76.59000  3.50 23.55000
## 186 46.51000       14.21   75.93000  3.58 24.07000
## 187 50.47000       18.67   77.59000  3.37 22.76000
## 188 51.38000       19.78   77.97000  3.32 24.83871
## 189 43.55000       11.28   74.70000  3.74 25.05000
## 190 44.01000       11.72   74.89000  3.72 24.89000
## 191 54.17000       23.38   79.14000  3.17 21.54000
## 193 43.01000       10.78   74.47000  3.77 25.22000
## 194 38.96000        7.40   72.78000  3.99 26.56000
## 195 55.94000       25.82   79.88000  3.07 20.95000
## 196 46.85000       14.57   76.08000  3.57 23.96000
## 197 37.95000        6.66   72.35000  4.05 26.90000
## 200 48.03000       15.85   76.57000  3.50 23.57000
## 201 47.11000       38.60   79.39516  3.50 23.87000
## 202 47.64000       15.42   76.41000  3.52 23.70000
## 203 44.08000       11.78   74.92000  3.72 24.87000
## 204 56.35000       13.34   78.64000  3.31 19.90000
## 205 51.58000        8.94   76.65000  3.57 21.48000
## 206 40.50000        2.10   72.02000  4.17 25.14000
## 208 67.73000       27.40   83.40000  2.69 16.15000
## 209 50.43000        8.01   76.17000  3.63 21.86000
## 211 54.67000       11.69   77.94000  3.40 20.46000
## 213 41.52000        2.53   72.44000  4.11 24.80000
## 214 48.20000        6.35   75.24000  3.75 22.60000
## 215 48.61000        6.64   75.41000  3.73 22.46000
## 216 49.93000        7.62   75.96000  3.66 22.03000
## 217 49.22000        7.08   75.66000  3.70 19.38710
## 218 45.59000        4.65   74.14000  3.89 23.46000
## 219 42.18000        2.83   72.72000  4.08 24.59000
## 220 56.43000       13.42   78.68000  3.31 19.88000
## 222 42.00000        2.75   72.64000  4.09 24.65000
## 223 40.50000        2.10   72.02000  4.17 25.14000
## 224 56.48000       13.47   78.70000  3.30 19.86000
## 225 47.53000        5.89   74.96000  3.79 22.82000
## 226 37.54000        1.08   70.78000  4.33 26.12000
## 229 50.50000        8.06   76.20000  3.63 21.84000
## 230 48.40000       17.49   78.45774  4.80 22.53000
## 231 43.46000        3.46   73.25000  4.01 24.16000
## 232 44.52000        4.03   73.70000  3.95 23.81000
## 233 51.14000        3.68   73.65000  3.51 21.21000
## 234 44.86000        1.09   71.03000  3.85 23.28000
## 235 39.24000        0.07   68.68000  4.16 25.14000
## 237 60.00000        9.91   77.36000  3.03 18.28000
## 238 42.91000        0.60   70.21000  3.96 23.92000
## 240 48.10000        2.23   72.38000  3.68 22.21000
## 242 41.36000        0.31   69.56000  4.04 24.44000
## 243 44.19000        0.90   70.75000  3.89 23.50000
## 244 42.80000        0.57   70.17000  3.96 23.96000
## 245 42.07000        0.43   69.86000  4.00 24.20000
## 246 43.47000        0.72   70.45000  3.93 20.44828
## 247 42.40000        0.49   70.00000  3.99 24.09000
## 248 38.36000        0.02   68.31000  4.20 25.43000
## 249 46.79000        1.72   71.84000  3.75 22.64000
## 251 35.40000        0.07   67.07000  4.36 26.41000
## 252 37.22000        0.00   67.83000  4.27 25.80000
## 253 53.40000        4.98   74.60000  3.39 20.46000
## 254 41.96000        0.41   69.82000  4.01 24.24000
## 255 37.02000        0.00   67.75000  4.28 25.87000
## 258 44.12000        0.89   70.72000  3.89 23.52000
## 259 46.53000        2.17   73.03300  4.30 22.73000
## 260 44.43000        0.97   70.85000  3.88 23.42000
## 261 42.35000        0.48   69.98000  3.99 24.11000
## 262 52.99000        9.57   73.18000  3.69 21.90000
## 263 58.55000       14.95   75.50000  3.38 20.07000
## 264 41.38000        2.19   68.32000  4.32 25.74000
## 266 62.27000       19.22   77.06000  3.18 18.84000
## 267 48.92000        6.39   71.48000  3.91 23.25000
## 269 52.56000        9.20   73.00000  3.71 22.04000
## 271 46.94000        5.07   70.65000  4.01 23.90000
## 272 50.34000        7.43   72.07000  3.83 22.78000
## 273 44.33000        3.57   69.56000  4.16 24.76000
## 274 48.69000        6.23   71.38000  3.92 23.32000
## 275 48.90000        6.38   71.47000  3.91 21.38710
## 276 45.97000        4.49   70.24000  4.07 24.22000
## 277 40.61000        1.88   68.00000  4.36 25.99000
## 278 53.80000       10.28   73.52000  3.64 21.64000
## 280 41.85000        2.39   68.52000  4.29 25.58000
## 281 40.48000        1.84   67.95000  4.36 26.04000
## 282 59.24000       15.70   75.79000  3.35 19.84000
## 283 46.57000        4.85   70.49000  4.03 24.02000
## 284 40.12000        1.70   67.80000  4.38 26.15000
## 287 47.90000        5.69   71.05000  3.96 23.58000
## 288 51.19000       34.50   71.84032  4.30 22.50000
## 289 48.43000        6.05   71.27000  3.93 23.41000
## 290 47.35000        5.33   70.82000  3.99 23.77000
## 291 42.65000      188.29   78.17000  3.08 23.90000
## 292 53.91000      233.69   82.88000  2.47 20.18000
## 293 36.90000      167.00   75.77000  3.39 25.80000
## 295 47.24000      206.21   80.09000  2.83 22.39000
## 296 40.11000      178.73   77.11000  3.21 24.74000
## 298 45.07000      197.64   79.18000  2.94 23.10000
## 300 40.48000      180.11   77.26000  3.19 24.62000
## 301 43.06000      189.86   78.34000  3.05 23.77000
## 302 36.82000      166.71   75.73000  3.39 25.83000
## 303 38.72000      173.60   76.53000  3.29 25.20000
## 304 42.66000      188.33   78.17000  3.08 29.62069
## 305 39.48000      176.40   76.84000  3.25 24.95000
## 306 30.90000      146.14   73.26000  3.71 27.79000
## 307 44.30000      194.64   78.86000  2.99 23.36000
## 309 35.58000      162.29   75.21000  3.46 26.24000
## 310 34.72000      159.26   74.85000  3.51 26.52000
## 311 51.31000      222.78   81.79000  2.61 21.04000
## 312 38.91000      174.30   76.61000  3.28 25.14000
## 313 35.04000      160.38   74.99000  3.49 26.42000
## 316 40.16000      178.92   77.13000  3.21 24.73000
## 317 41.58000      166.61   76.67733  3.40 24.26000
## 318 37.19000      168.04   75.89000  3.37 25.71000
## 319 41.50000      183.93   77.69000  3.14 24.28000
## 320 40.98000       89.25   77.99000  2.85  7.37000
## 321 45.70000      102.08   79.96000  2.59  5.81000
## 322 32.53000       68.44   74.46000  3.31 10.16000
## 324 39.30000       84.90   77.29000  2.94  7.92000
## 325 36.81000       78.64   76.25000  3.07  8.75000
## 327 42.16000       92.38   78.48000  2.78  6.98000
## 329 35.83000       76.24   75.84000  3.13  9.07000
## 330 38.32000       82.40   76.88000  2.99  8.25000
## 331 33.19000       69.97   74.73000  3.27  9.94000
## 332 42.57000       93.48   78.65000  2.76  6.84000
## 333 38.29000       82.33   76.87000  2.99 29.96774
## 334 33.86000       71.53   75.01000  3.23  9.72000
## 335 27.71000       57.81   72.44000  3.57 11.75000
## 336 41.84000       91.53   78.35000  2.80  7.08000
## 338 31.42000       65.91   73.99000  3.37 10.53000
## 339 30.78000       64.48   73.72000  3.40 10.74000
## 340 46.12000      103.26   80.14000  2.57  5.67000
## 341 35.08000       74.43   75.52000  3.17  9.32000
## 342 30.89000       64.72   73.77000  3.40 10.70000
## 345 35.72000       75.97   75.79000  3.13  9.11000
## 346 36.48000      128.27   77.53710  3.40  8.86000
## 347 30.72000       64.34   73.70000  3.40 10.76000
## 348 37.92000       81.40   76.71000  3.01  8.38000
## 
## Slot "sp":
## SpatialPoints:
##       coords.x1 coords.x2
##  [1,]  701078.1   9317300
##  [2,]  702363.8   9315200
##  [3,]  703520.7   9314565
##  [4,]  702072.5   9315722
##  [5,]  713562.7   9314648
##  [6,]  708202.7   9297513
##  [7,]  708365.4   9311897
##  [8,]  705238.6   9305209
##  [9,]  711732.0   9317081
## [10,]  711496.6   9321276
## [11,]  716938.0   9325901
## [12,]  712042.2   9316295
## [13,]  697455.1   9322826
## [14,]  706048.1   9320698
## [15,]  700097.1   9312272
## [16,]  702312.6   9304236
## [17,]  696675.0   9302755
## [18,]  694931.3   9305603
## [19,]  703922.7   9302476
## [20,]  693092.9   9317830
## [21,]  697690.2   9314998
## [22,]  688138.4   9317314
## [23,]  698100.5   9317878
## Coordinate Reference System (CRS) arguments: +proj=utm +zone=48 +south
## +datum=WGS84 +units=m +no_defs 
## 
## Slot "time":
##                     timeIndex
## 2023-01-01 12:00:00         1
## 2023-02-01 12:00:00         2
## 2023-03-01 12:00:00         3
## 2023-04-01 12:00:00         4
## 2023-05-01 12:00:00         5
## 2023-06-01 12:00:00         6
## 2023-07-01 12:00:00         7
## 2023-08-01 12:00:00         8
## 2023-09-01 12:00:00         9
## 2023-10-01 12:00:00        10
## 2023-11-01 12:00:00        11
## 2023-12-01 12:00:00        12
## 
## Slot "endTime":
##  [1] "2023-02-01 12:00:00 GMT" "2023-03-01 12:00:00 GMT"
##  [3] "2023-04-01 12:00:00 GMT" "2023-05-01 12:00:00 GMT"
##  [5] "2023-06-01 12:00:00 GMT" "2023-07-01 12:00:00 GMT"
##  [7] "2023-08-01 12:00:00 GMT" "2023-09-01 12:00:00 GMT"
##  [9] "2023-10-01 12:00:00 GMT" "2023-11-01 12:00:00 GMT"
## [11] "2023-12-01 12:00:00 GMT" "2023-12-31 12:00:00 GMT"

# STFDF DATA TEST
test_data = combined_data[combined_data$sp.ID %in% test_samples, ]
coords_test = unique(test_data[, c("coords.x1","coords.x2")])
pts_test = coordinates(coords_test)
pts_test = SpatialPoints(pts_test,
                         CRS("+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"))
utm48s = "+proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs"
pts_test.sfc = st_transform(st_as_sfc(pts_test), utm48s)
pts_test = as(pts_test.sfc, "Spatial")
stfdf.test = STFDF(sp = pts_test, time = time, data = test_data[,7:11])
stfdf.test

## An object of class "STFDF"
## Slot "data":
##      PM25 Presipitasi Kelembaban Angin   NO2
## 4   38.08      219.27      85.86  4.17  5.13
## 7   18.34      145.54      77.61  5.24 11.65
## 9   27.43      177.62      81.41  4.75  8.65
## 18  22.42      159.54      79.31  5.02 10.30
## 24  36.73      213.74      85.30  4.25  5.58
## 25  30.50      189.17      82.69  4.58  7.63
## 33  28.26      549.53      87.18  5.48 11.92
## 36  15.72      470.85      81.94  6.15 16.07
## 38  24.23      523.59      85.50  5.69 13.25
## 47  18.64      488.63      83.16  6.00 15.10
## 53  29.68      558.83      87.78  5.40 11.45
## 54  29.52      557.78      87.71  5.41 11.51
## 62  46.60      317.48      88.97  3.17 18.02
## 65  26.97      227.70      80.76  4.24 24.50
## 67  29.17      237.02      81.68  4.12 23.78
## 76  29.54      238.61      81.83  4.10 23.65
## 82  46.47      316.84      88.91  3.18 18.06
## 83  48.57      327.32      89.79  3.07 17.37
## 91  41.02      215.53      84.48  3.01 12.98
## 94  24.76      154.28      77.68  3.89 18.35
## 96  26.41      160.03      78.37  3.80 17.81
## 105 30.35      174.18      80.02  3.59 16.51
## 111 35.21      192.47      82.05  3.33 14.90
## 112 34.84      191.05      81.90  3.35 15.02
## 120 64.90      199.50      87.96  2.21 17.72
## 123 40.47      115.11      77.74  3.53 25.79
## 125 37.37      106.05      76.44  3.70 26.81
## 134 50.78      147.91      82.05  2.97 22.38
## 140 58.24      174.21      85.17  2.57 19.92
## 141 41.48      118.14      78.16  3.48 25.45
## 149 61.65      189.58      87.95  2.14 15.26
## 152 39.84      115.29      78.83  3.32 22.47
## 154 36.99      106.94      77.64  3.48 23.41
## 163 49.96      147.47      83.06  2.77 19.13
## 169 53.64      160.16      84.60  2.57 17.91
## 170 58.95      179.38      86.82  2.29 16.16
## 178 64.58       39.46      83.49  2.61 18.10
## 181 42.36       10.20      74.20  3.81 25.44
## 183 41.01        9.03      73.63  3.88 25.89
## 192 53.05       21.90      78.67  3.23 21.91
## 198 56.32       26.36      80.04  3.05 20.83
## 199 64.97       40.15      83.65  2.58 17.97
## 207 68.83       29.02      83.86  2.63 15.78
## 210 42.95        3.20      73.04  4.04 24.33
## 212 40.98        2.30      72.22  4.14 24.98
## 221 52.30        9.54      76.95  3.53 21.24
## 227 58.65       15.78      79.61  3.19 19.15
## 228 67.18       26.60      83.17  2.72 16.33
## 236 62.13       11.87      78.25  2.92 17.58
## 239 41.46        0.33      69.61  4.04 24.40
## 241 39.94        0.13      68.97  4.12 24.91
## 250 48.13        2.25      72.40  3.67 22.20
## 256 53.15        4.83      74.49  3.40 20.54
## 257 62.53       12.25      78.42  2.89 17.44
## 265 65.71       23.64      78.50  3.00 17.70
## 268 41.12        2.08      68.21  4.33 25.82
## 270 44.18        3.49      69.49  4.16 24.81
## 279 54.72       11.12      73.90  3.59 21.33
## 285 61.22       17.96      76.62  3.24 19.18
## 286 68.15       27.05      79.52  2.86 16.89
## 294 54.38      235.70      83.07  2.44 20.03
## 297 33.43      154.77      74.31  3.58 26.95
## 299 38.03      171.08      76.24  3.33 25.43
## 308 44.16      194.10      78.80  2.99 23.40
## 314 53.41      231.57      82.67  2.49 20.35
## 315 54.37      235.65      83.07  2.44 20.03
## 323 51.49      118.99      82.38  2.28  3.90
## 326 30.51       63.87      73.61  3.42 10.83
## 328 34.24       72.43      75.17  3.21  9.60
## 337 40.02       86.75      77.59  2.90  7.69
## 343 45.79      102.33      80.00  2.59  5.78
## 344 47.62      107.54      80.77  2.49  5.18
## 
## Slot "sp":
## SpatialPoints:
##      coords.x1 coords.x2
## [1,]  702838.7   9315468
## [2,]  706282.9   9308147
## [3,]  712173.4   9301736
## [4,]  699318.8   9306375
## [5,]  690991.0   9322319
## [6,]  699824.2   9319688
## Coordinate Reference System (CRS) arguments: +proj=utm +zone=48 +south
## +datum=WGS84 +units=m +no_defs 
## 
## Slot "time":
##                     timeIndex
## 2023-01-01 12:00:00         1
## 2023-02-01 12:00:00         2
## 2023-03-01 12:00:00         3
## 2023-04-01 12:00:00         4
## 2023-05-01 12:00:00         5
## 2023-06-01 12:00:00         6
## 2023-07-01 12:00:00         7
## 2023-08-01 12:00:00         8
## 2023-09-01 12:00:00         9
## 2023-10-01 12:00:00        10
## 2023-11-01 12:00:00        11
## 2023-12-01 12:00:00        12
## 
## Slot "endTime":
##  [1] "2023-02-01 12:00:00 GMT" "2023-03-01 12:00:00 GMT"
##  [3] "2023-04-01 12:00:00 GMT" "2023-05-01 12:00:00 GMT"
##  [5] "2023-06-01 12:00:00 GMT" "2023-07-01 12:00:00 GMT"
##  [7] "2023-08-01 12:00:00 GMT" "2023-09-01 12:00:00 GMT"
##  [9] "2023-10-01 12:00:00 GMT" "2023-11-01 12:00:00 GMT"
## [11] "2023-12-01 12:00:00 GMT" "2023-12-31 12:00:00 GMT"

Empirical Semivariogram

# membuat object gstat untuk variogram
g <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("Presipitasi", Presipitasi ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("Kelembaban", Kelembaban ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("Angin", Angin ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("NO2", NO2 ~ 1, data = stfdf.train, fill.cross = TRUE)

# =========================================Menghitung cross-variogram empiris
var <- variogramST(g, data = stfdf.train, 
                   assumeRegular = TRUE, pseudo = 1)

#write.xlsx(var, "D:/TESIS/Hasil/excels/Semivariogram empiris.xlsx")

# Ekstrak informasi dari variogramST object
gamma_values <- var$gamma  # Nilai variogram (semivariance)
distances <- var$spacelag  # Jarak (spatial lag)
time_lags <- var$timelag   # Time lag (temporal)

# Buat plot interaktif 3D
plot_ly(x = ~distances, y = ~time_lags, z = ~gamma_values, 
        intensity = ~gamma_values, type = 'mesh3d', 
        colors = colorRamp(c("darkblue", "blue", "red", "orange", "yellow"))) %>%
  layout(scene = list(
    xaxis = list(title = "Jarak spasial"),
    yaxis = list(title = "Lag Waktu (hari)"),
    zaxis = list(title = "Gamma")
  )) %>%
  colorbar(title = "Gamma")

Semivariogram empiris dihasilkan menggunakan data pelatihan yang terdiri dari 23 titik lokasi selama periode 12 bulan. Komponen spasial dari semivariogram empiris menunjukkan variasi minimal dalam nilai gamma seiring dengan bertambahnya jarak spasial, yang mengindikasikan variabilitas spasial yang terbatas dalam rentang pengamatan. Sebaliknya, nilai gamma menunjukkan peningkatan yang jelas seiring dengan bertambahnya lag waktu, mencerminkan variabilitas yang lebih besar seiring waktu. Menariknya, penurunan nilai gamma diamati pada lag waktu terbesar, yang mungkin menunjukkan efek smoothing atau korelasi temporal yang berkurang pada interval yang lebih panjang.

Semivariogram empiris kemudian dimodelkan menggunakan berbagai semivariogram teoretis. Dengan menggabungkan tiga model semivariogram marginal (spherical, exponential, dan Gaussian) untuk komponen spasial, temporal, dan gabungan (spasial), total 9 model separable, 9 model product-sum, 3 model metric, 27 model sum-metric, dan 27 model simple sum-metric dievaluasi. Gambar menggambarkan bentuk semivariogram spasiotemporal yang dihasilkan dari kombinasi terbaik di setiap kategori model teoretis. Proses pemodelan ini bertujuan untuk menentukan semivariogram teoretis mana yang paling akurat menangkap pola variabilitas spasiotemporal yang diamati dalam semivariogram empiris.

head(var, 6)

##   np      dist    gamma   id timelag spacelag   avgDist
## 1  0        NA       NA lag0  0 days    0.000    0.0000
## 2 24  721.5611 2.119802 lag0  0 days  449.318  721.5611
## 3 12 1319.6596 8.503237 lag0  0 days 1347.954 1319.6596
## 4 60 2161.6545 4.327761 lag0  0 days 2246.590 2161.6545
## 5 72 3074.1309 2.719611 lag0  0 days 3145.226 3074.1309
## 6 96 3980.1009 6.023574 lag0  0 days 4043.862 3980.1009

tail(var, 6)

##     np      dist     gamma    id  timelag  spacelag   avgDist
## 187 26  8604.812  2.860425 lag11 341 days  8537.042  8604.812
## 188 28  9515.440  3.537421 lag11 341 days  9435.678  9515.440
## 189 34 10284.410 17.669637 lag11 341 days 10334.314 10284.410
## 190 26 11077.128  2.925169 lag11 341 days 11232.950 11077.128
## 191 22 12132.165 30.303786 lag11 341 days 12131.586 12132.165
## 192 30 12892.806  2.563838 lag11 341 days 13030.222 12892.806

# setting lower and upper bounds
pars.l <- c(sill.s = 13.75739, range.s = 0.1*30762.69, nugget.s = 0, 
            sill.t = 43.72425, range.t = 31, nugget.t = 0, 
            sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0, 
            anis = 0)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745, 
            sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959, 
            sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855, 
            anis = 928) 

# separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     separable <- vgmST("separable", space = vgm(psill = 13.75739, model = sm, range = 0.1*30762.69, nugget = 0),
#                        time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 43.72425)
#     separable_Vgm <- fit.StVariogram(var, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
#     MSE = c(MSE, attr(separable_Vgm, "MSE"))
#   }
# }
# separable = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]
# write.xlsx(separable, "D:/TESIS/Hasil/excels/MSE Separable.xlsx")

# product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                           time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 0.5) 
#     # Menggunakan tryCatch untuk menangani error
#     fit_result <- tryCatch({
#       fit.StVariogram(var, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#     }, error = function(e) {
#       return(NULL)  # Jika ada error, kembalikan NULL
#     })

#     # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#     if (is.null(fit_result)) {
#       MSE = c(MSE, 0)
#     } else {
#       MSE = c(MSE, attr(fit_result, "MSE"))
#     }
#   }
# }
# prodsum = data.frame(Spatial = model_space, 
#                      Temporal = model_temporal,
#                      MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]
# write.xlsx(prodsum, "D:/TESIS/Hasil/excels/MSE Product-Sum.xlsx")

# Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
#     model_joint = c(model_joint, j)
#     metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = j, range = 0.5*30762.69, nugget = 0), stAni=1) 
#     metric_Vgm <- fit.StVariogram(var, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
#     MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
#                     MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]
# write.xlsx(metric, "D:/TESIS/Hasil/excels/MSE Metric.xlsx")

# Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                          time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                          joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(var, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })

#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# summetric = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        Joint = model_joint,
#                        MSE = as.numeric(MSE))
# summetric
# write.xlsx(summetric, "D:/TESIS/Hasil/excels/MSE Sum Metric.xlsx")

# Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                                time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                                joint = vgm(psill = 136.6855, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(var, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })

#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# ssmetric = data.frame(Spatial = model_space,
#                       Temporal = model_temporal,
#                       Joint = model_joint,
#                       MSE = as.numeric(MSE))
# ssmetric
# write.xlsx(ssmetric, "D:/TESIS/Hasil/excels/MSE Simple-Sum Metric.xlsx")

# separable = read.xlsx("D:/TESIS/Hasil/excels/MSE Separable.xlsx")
# prodsum = read.xlsx("D:/TESIS/Hasil/excels/MSE Product-Sum.xlsx")
# metric = read.xlsx("D:/TESIS/Hasil/excels/MSE Metric.xlsx")
# summetric = read.xlsx("D:/TESIS/Hasil/excels/MSE Sum Metric.xlsx")
# ssmetric = read.xlsx("D:/TESIS/Hasil/excels/MSE Simple-Sum Metric.xlsx")

# separable[separable$MSE==min(separable$MSE),]
# prodsum=prodsum[prodsum$MSE>0,]
# prodsum[prodsum$MSE==min(prodsum$MSE, na.rm = T),]
# metric[metric$MSE==min(metric$MSE, na.rm = T),]
# summetric=summetric[summetric$MSE>0,]
# summetric[summetric$MSE==min(summetric$MSE, na.rm = T),]
# ssmetric=ssmetric[ssmetric$MSE>0,]
# ssmetric[ssmetric$MSE==min(ssmetric$MSE, na.rm = T),]

# ==============================================================MODEL TERBAIK
# separable
separable <- vgmST("separable", space = vgm(psill = 13.75739, model = 'Exp', range = 0.5*30762.69, nugget = 0),
                   time = vgm(psill = 109.7589, model = 'Exp', range = 217, nugget = 0), sill = 43.72425)
separable_Vgm <- fit.StVariogram(var, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
param_sep = extractPar(separable_Vgm)

# product-sum
prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
                      time = vgm(psill = 109.7589, model = 'Gau', range = 217, nugget = 0), k = 10) 
prodSumModel_Vgm <- fit.StVariogram(var, prodSumModel, method = "L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
param_prodsum = extractPar(prodSumModel_Vgm)

# Metric
metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = 'Gau', range = 0.5*30762.69, nugget = 0), stAni=1.5) 
metric_Vgm <- fit.StVariogram(var, metric, method="L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
param_metric = extractPar(metric_Vgm)

# Sum-Metric
sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
                   time = vgm(psill = 109.7589, model = 'Exp', range = 217, nugget = 0), 
                   joint = vgm(psill = 136.6855, model = 'Sph', range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
sumMetric_Vgm <- fit.StVariogram(var, sumMetric, method="L-BFGS-B", tunit="days", fit.method = 6, lower = pars.l, upper = pars.u)
param_summetric = extractPar(sumMetric_Vgm)

# Simple Sum-Metric
SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
                         time = vgm(psill = 109.7589, model = 'Gau', range = 217, nugget = 0), 
                         joint = vgm(psill = 136.6855, model = 'Exp', range = sqrt((31000^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=0.5) 
SimplesumMetric_Vgm <- fit.StVariogram(var, SimplesumMetric, method = "L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
param_ssmetric = extractPar(SimplesumMetric_Vgm)

# =========================================================PLOT SEMIVARIOGRAM
sampel_plot = plot(var, wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
sep_plot = plot(var,separable_Vgm, wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
ps_plot = plot(var,prodSumModel_Vgm,wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
metric_plot = plot(var, metric_Vgm, wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
summet_plot = plot(var,sumMetric_Vgm,wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)
ssmet_plot = plot(var,SimplesumMetric_Vgm,wireframe=T, xlab = NULL, ylab = NULL, zlab = NULL)

sampel_plot = as.ggplot(sampel_plot)
sep_plot <- as.ggplot(sep_plot)
ps_plot <- as.ggplot(ps_plot)
metric_plot <- as.ggplot(metric_plot)
summet_plot <- as.ggplot(summet_plot)
ssmet_plot <- as.ggplot(ssmet_plot)

# Combine the plots using wrap_plots
semivar_plot <- wrap_plots(
  list(sampel_plot, sep_plot, ps_plot, metric_plot, summet_plot, ssmet_plot), 
  ncol = 3) +
  plot_layout(guides = "auto") +
  plot_annotation(title = "Perbandingan Semivariogram Empiris dengan Semivariogram Teoretis")

# Print the combined plot
print(semivar_plot)

Semivariogram teoretis dengan MSE terendah menunjukkan kemampuannya untuk menjelaskan pola variabilitas yang diamati dalam semivariogram empiris. Oleh karena itu, semivariogram teoretis dengan MSE terendah akan digunakan untuk menghasilkan matriks gamma dalam persamaan prediksi cokriging. Model teoretis dengan MSE terendah adalah model sum-metric dengan nilai MSE sebesar 292,152. Selain itu, juga terlihat bahwa semivariogram sum-metric menunjukkan pola yang paling mirip dengan semivariogram empiris. Dengan demikian, model sum-metric akan digunakan untuk menghasilkan matriks gamma \(\Gamma_{11}\) untuk perhitungan bobot cokriging selanjutnya.

Cross-Variogram

Proses fitting juga dilakukan terhadap cross-variogram antara variabel \(PM_{2.5}\) dengan variabel prediktor menggunakan langkah yang sama seperti pada fitting semivariogram. Fitting dilakukan menggunakan model teoretis product, product-sum, metric, sum-metric, dan simple sum-metric. Proses ini dilakukan untuk mendapatkan estimasi nilai parameter cross-variogram dengan meminimalkan selisih antara nilai empiris dengan teoretis menggunakan metode OLS. Proses ini akan menghasilkan fungsi kontinu yang akan digunakan dalam pembangkitan nilai cross-variogram.

Cross-Variogram antara \(PM_{2.5}\) and \(NO_{2}\)

# =========================================PM25~NO2
g_2 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("NO2", NO2 ~ 1, data = stfdf.train, fill.cross = TRUE)

# Menghitung cross-variogram empiris
cross_var_2 <- variogramST(g_2, data = stfdf.train, 
                           assumeRegular = TRUE, pseudo = 1)

pars.l <- c(sill.s = 13.75739, range.s = 0.1*30762.69, nugget.s = 0, 
            sill.t = 43.72425, range.t = 31, nugget.t = 0, 
            sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0, 
            anis = 10)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745, 
            sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959, 
            sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855, 
            anis = 928)
# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     separable <- vgmST("separable", space = vgm(psill = 13.75739, model = sm, range = 0.1*30762.69, nugget = 0),
#                        time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 43.72425)
#     separable_Vgm <- fit.StVariogram(cross_var_2, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
#     MSE = c(MSE, attr(separable_Vgm, "MSE"))
#   }
# }
# separable = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]

# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                           time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 0.5) 
#     # Menggunakan tryCatch untuk menangani error
#     fit_result <- tryCatch({
#       fit.StVariogram(cross_var_2, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#     }, error = function(e) {
#       return(NULL)  # Jika ada error, kembalikan NULL
#     })
#     
#     # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#     if (is.null(fit_result)) {
#       MSE = c(MSE, 0)
#     } else {
#       MSE = c(MSE, attr(fit_result, "MSE"))
#     }
#   }
# }
# prodsum = data.frame(Spatial = model_space, 
#                      Temporal = model_temporal,
#                      MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]

# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
#   model_joint = c(model_joint, j)
#   metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = j, range = 0.5*30762.69, nugget = 0), stAni=1) 
#   metric_Vgm <- fit.StVariogram(cross_var_2, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
#   MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
#                     MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]

# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                          time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                          joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_2, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# summetric = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        Joint = model_joint,
#                        MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]

# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                                time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                                joint = vgm(psill = 136.6855, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_2, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# ssmetric = data.frame(Spatial = model_space,
#                       Temporal = model_temporal,
#                       Joint = model_joint,
#                       MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]
# separable$Model = c(rep("Product",9))
# separable$Joint = c(rep(" ",9))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))

# pm_no = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_no, "D:/TESIS/Hasil/Cross-Variogram PM25 dan NO2.xlsx")

# Sum-Metric
sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
                   time = vgm(psill = 109.7589, model = 'Sph', range = 217, nugget = 0), 
                   joint = vgm(psill = 136.6855, model = 'Exp', range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=1.5)
sumMetric_Vgm <- fit.StVariogram(cross_var_2, sumMetric, method="L-BFGS-B", tunit="days", fit.method = 6, lower = pars.l, upper = pars.u)
attr(sumMetric_Vgm, "MSE")

## [1] 288.9448

cv2_summetric = extractPar(sumMetric_Vgm);cv2_summetric

##       sill.s      range.s     nugget.s       sill.t      range.t     nugget.t 
## 1.375739e+01 3.076269e+04 0.000000e+00 4.372425e+01 2.809691e+02 2.827152e+00 
##      sill.st     range.st    nugget.st         anis 
## 1.366855e+01 2.195445e+04 5.637747e-01 1.000000e+01

plot(cross_var_2, list(sumMetric_Vgm), all=T, wireframe=T)

Pada fitting ini, fungsi semivariogram marginal spasial menggunakan model Gaussian, sementara fungsi marginal temporal menggunakan model spherical, dan komponen joint menggunakan fungsi eksponensial. Gambar tersebut menunjukkan perbandingan cross-variogram empiris variabel konsentrasi \(PM_{2.5}\) dan konsentrasi \(NO_2\) dengan cross-variogram teoretis model sum-metric.

Cross-Variogram antara \(PM_{2.5}\) and Curah Hujan

# =========================================PM25~Presipitasi
g_3 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("Presipitasi", Presipitasi ~ 1, data = stfdf.train, fill.cross = TRUE)

# Menghitung cross-variogram empiris
cross_var_3 <- variogramST(g_3, data = stfdf.train, 
                           assumeRegular = TRUE, pseudo = 1)

pars.l <- c(sill.s = 914.197, range.s = 0.1*30762.69, nugget.s = 0, 
            sill.t = 43.72425, range.t = 31, nugget.t = 0, 
            sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0, 
            anis = 10)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745, 
            sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959, 
            sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855, 
            anis = 928)

# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     separable <- vgmST("separable", space = vgm(psill = 100, model = sm, range = 0.1*30762.69, nugget = 0),
#                        time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 50)
#     separable_Vgm <- fit.StVariogram(cross_var_3, separable, fit.method = 6, method="L-BFGS-B")
#     MSE = c(MSE, attr(separable_Vgm, "MSE"))
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#   }
# }
# separable = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]

# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                           time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 10) 
#     # Menggunakan tryCatch untuk menangani error
#     fit_result <- tryCatch({
#       fit.StVariogram(cross_var_3, prodSumModel, fit.method = 6, method = "L-BFGS-B")
#     }, error = function(e) {
#       return(NULL)  # Jika ada error, kembalikan NULL
#     })
#     
#     # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#     if (is.null(fit_result)) {
#       MSE = c(MSE, 0)
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#     } else {
#       MSE = c(MSE, attr(fit_result, "MSE"))
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#     }
#   }
# }
# prodsum = data.frame(Spatial = model_space, 
#                      Temporal = model_temporal,
#                      MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]

# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
#   model_joint = c(model_joint, j)
#   metric <- vgmST("metric", joint = vgm(psill = 45000, model = j, range = 0.5*30762.69, nugget = 0), stAni=1.5) 
#   metric_Vgm <- fit.StVariogram(cross_var_3, metric, method="L-BFGS-B")
#   MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
#                     MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]

# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                          time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                          joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=1.5)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_3, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#         model_space = c(model_space, sm)
#         model_temporal = c(model_temporal, tm)
#         model_joint = c(model_joint, j)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#         model_space = c(model_space, sm)
#         model_temporal = c(model_temporal, tm)
#         model_joint = c(model_joint, j)
#       }
#     }
#   }
# }
# summetric = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        Joint = model_joint,
#                        MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]

# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                                time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                                joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=0.5)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_3, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#         model_space = c(model_space, sm)
#         model_temporal = c(model_temporal, tm)
#         model_joint = c(model_joint, j)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#         model_space = c(model_space, sm)
#         model_temporal = c(model_temporal, tm)
#         model_joint = c(model_joint, j)
#       }
#     }
#   }
# }
# ssmetric = data.frame(Spatial = model_space,
#                       Temporal = model_temporal,
#                       Joint = model_joint,
#                       MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]

# separable$Model = c(rep("Product", 6))
# separable$Joint = c(rep(" ", 6))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))

# pm_ch = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_ch, "D:/TESIS/Hasil/Cross-Variogram PM25 dan Curah Hujan.xlsx")

# Simple Sum-Metric
SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = 'Gau', range = 0.5*30762.69, nugget = 0),
                         time = vgm(psill = 109.7589, model = 'Sph', range = 217, nugget = 0), 
                         joint = vgm(psill = 136.6855, model = 'Gau', range = sqrt((31000^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=0.5) 
SimplesumMetric_Vgm <- fit.StVariogram(cross_var_3, SimplesumMetric, method = "L-BFGS-B", fit.method = 6, lower = pars.l)
attr(SimplesumMetric_Vgm, "MSE")

## [1] 83946136

cv3_ssmetric = extractPar(SimplesumMetric_Vgm);cv3_ssmetric

##      sill.s     range.s      sill.t     range.t     sill.st    range.st 
##    914.1970 659784.4230  26804.7595    238.9169   1718.8804  49970.7851 
##      nugget        anis 
##    302.3135   9728.0167

plot(cross_var_3, list(SimplesumMetric_Vgm), all=T, wireframe=T)

Dalam fitting cross-variogram antara variabel konsentrasi \(PM_{2.5}\) dan curah hujan, model simple sum-metric menghasilkan nilai MSE paling rendah dibandingkan model lainnya. Model tersebut menggunakan fungsi Gaussian pada komponen spasial, fungsi spherical pada temporal, dan fungsi Gaussian pada komponen joint.

Cross-Variogram antara \(PM_{2.5}\) and Air Humidity$

# =========================================PM25~Kelembaban
g_4 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("Kelembaban", Kelembaban ~ 1, data = stfdf.train, fill.cross = TRUE)

# Menghitung cross-variogram empiris
cross_var_4 <- variogramST(g_4, data = stfdf.train, 
                           assumeRegular = TRUE, pseudo = 1)

pars.l <- c(sill.s = 13.75739, range.s = 0.1*30762.69, nugget.s = 0, 
            sill.t = 43.72425, range.t = 31, nugget.t = 0, 
            sill.st = 0.1*136.6855, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0, 
            anis = 10)
pars.u <- c(sill.s = 70.55745, range.s = 30762.69, nugget.s = 0.5*70.55745, 
            sill.t = 214.7959, range.t = 341, nugget.t = 0.5*214.7959, 
            sill.st = 136.6855, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855, 
            anis = 928)

# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     separable <- vgmST("separable", space = vgm(psill = 13.75739, model = sm, range = 0.1*30762.69, nugget = 0),
#                        time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), sill = 43.72425)
#     separable_Vgm <- fit.StVariogram(cross_var_4, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
#     MSE = c(MSE, attr(separable_Vgm, "MSE"))
#   }
# }
# separable = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]

# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     prodSumModel <- vgmST("productSum", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                           time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), k = 0.5) 
#     # Menggunakan tryCatch untuk menangani error
#     fit_result <- tryCatch({
#       fit.StVariogram(cross_var_4, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#     }, error = function(e) {
#       return(NULL)  # Jika ada error, kembalikan NULL
#     })
#     
#     # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#     if (is.null(fit_result)) {
#       MSE = c(MSE, 0)
#     } else {
#       MSE = c(MSE, attr(fit_result, "MSE"))
#     }
#   }
# }
# prodsum = data.frame(Spatial = model_space, 
#                      Temporal = model_temporal,
#                      MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]

# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
#   model_joint = c(model_joint, j)
#   metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = j, range = 0.5*30762.69, nugget = 0), stAni=1) 
#   metric_Vgm <- fit.StVariogram(cross_var_4, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
#   MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
#                     MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]

# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       sumMetric <- vgmST("sumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                          time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                          joint = vgm(psill = 136.6855, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_4, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# summetric = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        Joint = model_joint,
#                        MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]

# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 13.75739, model = sm, range = 0.5*30762.69, nugget = 0),
#                                time = vgm(psill = 109.7589, model = tm, range = 217, nugget = 0), 
#                                joint = vgm(psill = 136.6855, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_4, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# ssmetric = data.frame(Spatial = model_space,
#                       Temporal = model_temporal,
#                       Joint = model_joint,
#                       MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]

# separable$Model = c(rep("Product",9))
# separable$Joint = c(rep(" ",9))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))

# pm_ku = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_ku, "D:/TESIS/Hasil/Cross-Variogram PM25 dan Kelembaban.xlsx")

# Metric
metric <- vgmST("metric", joint = vgm(psill = 136.6855, model = 'Gau', range = 0.5*30762.69, nugget = 0), stAni=1.5) 
metric_Vgm <- fit.StVariogram(cross_var_4, metric, method="L-BFGS-B", fit.method = 6, lower = pars.l, upper = pars.u)
attr(metric_Vgm, "MSE")

## [1] 78.54296

cv4_metric = extractPar(metric_Vgm);cv4_metric

##         sill        range       nugget         anis 
##    24.427639 22603.546735     6.308559   214.795900

plot(cross_var_4, list(metric_Vgm), all=T, wireframe=T, xlab = "Jarak Spasial", ylab = 'Lag Waktu (hari)', 
     zlab = 'Gamma')

Hasil fitting cross-variogram variabel konsentrasi \(PM_{2.5}\) dengan kelembaban menunjukkan bahwa model metric dengan fungsi Gaussian memberikan MSE terendah dibandingkan model lainnya

Cross-Variogram antara \(PM_{2.5}\) and Wind Speed$

# =========================================PM25~Angin
g_5 <- gstat(NULL, "PM25", PM25 ~ 1, data = stfdf.train, fill.cross = TRUE) %>%
  gstat("Angin", Angin ~ 1, data = stfdf.train, fill.cross = TRUE)

# Menghitung cross-variogram empiris
cross_var_5 <- variogramST(g_5, data = stfdf.train, 
                           assumeRegular = TRUE, pseudo = 1)

pars.l <- c(sill.s = 1, range.s = 0.1*30762.69, nugget.s = 0, 
            sill.t = 1, range.t = 31, nugget.t = 0, 
            sill.st = 1, range.st = sqrt((0.1*30762.69^2) + (31*0)^2), nugget.st = 0, 
            anis = 10)
pars.u <- c(sill.s = 3, range.s = 30762.69, nugget.s = 0.5*70.55745, 
            sill.t = 3, range.t = 341, nugget.t = 0.5*214.7959, 
            sill.st = 3, range.st = sqrt((30762.69^2) + (341*928)^2), nugget.st = 0.5*136.6855, 
            anis = 928)

# # separable
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     separable <- vgmST("separable", space = vgm(psill = 1, model = sm, range = 0.1*30762.69, nugget = 0),
#                        time = vgm(psill = 1, model = tm, range = 217, nugget = 0), sill = 1)
#     separable_Vgm <- fit.StVariogram(cross_var_5, separable, fit.method = 6, method="L-BFGS-B", upper = pars.u)
#     MSE = c(MSE, attr(separable_Vgm, "MSE"))
#   }
# }
# separable = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        MSE = as.numeric(MSE))
# separable[separable$MSE==min(separable$MSE),]

# # product-sum
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     model_space = c(model_space, sm)
#     model_temporal = c(model_temporal, tm)
#     prodSumModel <- vgmST("productSum", space = vgm(psill = 1, model = sm, range = 0.5*30762.69, nugget = 0),
#                           time = vgm(psill = 1, model = tm, range = 217, nugget = 0), k = 0.5) 
#     # Menggunakan tryCatch untuk menangani error
#     fit_result <- tryCatch({
#       fit.StVariogram(cross_var_5, prodSumModel, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#     }, error = function(e) {
#       return(NULL)  # Jika ada error, kembalikan NULL
#     })
#     
#     # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#     if (is.null(fit_result)) {
#       MSE = c(MSE, 0)
#     } else {
#       MSE = c(MSE, attr(fit_result, "MSE"))
#     }
#   }
# }
# prodsum = data.frame(Spatial = model_space, 
#                      Temporal = model_temporal,
#                      MSE = as.numeric(MSE))
# prodsum[prodsum$MSE==min(prodsum$MSE),]

# # Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_joint = c()
# MSE = c()
# for(j in daftar_model){
#   model_joint = c(model_joint, j)
#   metric <- vgmST("metric", joint = vgm(psill = 3, model = j, range = 0.5*30762.69, nugget = 0), stAni=1) 
#   metric_Vgm <- fit.StVariogram(cross_var_5, metric, method="L-BFGS-B", lower = pars.l, upper = pars.u)
#   MSE = c(MSE, attr(metric_Vgm, "MSE"))
# }
# metric = data.frame(Joint = model_joint,
#                     MSE = as.numeric(MSE))
# metric[metric$MSE==min(metric$MSE),]

# # Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       sumMetric <- vgmST("sumMetric", space = vgm(psill = 1, model = sm, range = 0.5*30762.69, nugget = 0),
#                          time = vgm(psill = 1, model = tm, range = 217, nugget = 0), 
#                          joint = vgm(psill = 3, model = j, range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_5, sumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# summetric = data.frame(Spatial = model_space, 
#                        Temporal = model_temporal,
#                        Joint = model_joint,
#                        MSE = as.numeric(MSE))
# summetric[summetric$MSE==min(summetric$MSE),]

# # Simple Sum-Metric
# daftar_model = c("Gau", "Sph", "Exp")
# model_space = c()
# model_temporal = c()
# model_joint = c()
# MSE = c()
# for(sm in daftar_model){
#   for(tm in daftar_model){
#     for(j in daftar_model){
#       model_space = c(model_space, sm)
#       model_temporal = c(model_temporal, tm)
#       model_joint = c(model_joint, j)
#       SimplesumMetric <- vgmST("simpleSumMetric", space = vgm(psill = 1, model = sm, range = 0.5*30762.69, nugget = 0),
#                                time = vgm(psill = 1, model = tm, range = 217, nugget = 0), 
#                                joint = vgm(psill = 3, model = j, range = sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), nugget = 0, stAni=500, k = 0.5)
#       # Menggunakan tryCatch untuk menangani error
#       fit_result <- tryCatch({
#         fit.StVariogram(cross_var_5, SimplesumMetric, fit.method = 6, method = "L-BFGS-B", lower = pars.l, upper = pars.u)
#       }, error = function(e) {
#         return(NULL)  # Jika ada error, kembalikan NULL
#       })
#       
#       # Jika hasil fitting ada error, tambahkan NA ke MSE, jika tidak, tambahkan MSE dari model
#       if (is.null(fit_result)) {
#         MSE = c(MSE, 0)
#       } else {
#         MSE = c(MSE, attr(fit_result, "MSE"))
#       }
#     }
#   }
# }
# ssmetric = data.frame(Spatial = model_space,
#                       Temporal = model_temporal,
#                       Joint = model_joint,
#                       MSE = as.numeric(MSE))
# ssmetric[ssmetric$MSE==min(ssmetric$MSE),]

# separable$Model = c(rep("Product",9))
# separable$Joint = c(rep(" ",9))
# prodsum$Model = c(rep("Product-Sum", 9))
# prodsum$Joint = c(rep(" ", 9))
# metric$Model = c(rep("Model", 3))
# metric$Spatial = c(rep(" ", 3))
# metric$Temporal = c(rep(" ", 3))
# summetric$Model = c(rep("Sum-Metric", 27))
# ssmetric$Model = c(rep("Simple Sum-Metric", 27))

# pm_angin = dplyr::bind_rows(separable, prodsum, metric, summetric, ssmetric)
# write.xlsx(pm_angin, "D:/TESIS/Hasil/Cross-Variogram PM25 dan Angin.xlsx")

# Sum-Metric
sumMetric <- vgmST("sumMetric", space = vgm(psill = 1, model = "Gau", range = 0.5*30762.69, nugget = 0),
                   time = vgm(psill = 1, model = "Gau", range = 217, nugget = 0), 
                   joint = vgm(psill = 3, model = "Gau", range = 0.5*sqrt((30762.69^2) + (341*1.5)^2), nugget = 0), stAni=500)
sumMetric_Vgm <- fit.StVariogram(cross_var_5, sumMetric, method="L-BFGS-B", tunit="days", fit.method = 6, lower = pars.l, upper = pars.u)
attr(sumMetric_Vgm, "MSE")

## [1] 0.1726448

cv5_summetric = extractPar(sumMetric_Vgm);cv5_summetric

##       sill.s      range.s     nugget.s       sill.t      range.t     nugget.t 
## 1.000000e+00 3.076269e+04 0.000000e+00 1.883009e+00 3.410000e+02 2.155019e-01 
##      sill.st     range.st    nugget.st         anis 
## 1.000000e+00 2.684337e+04 0.000000e+00 6.152708e+01

plot(cross_var_5, list(sumMetric_Vgm), all=T, wireframe=T, xlab = "Jarak Spasial", ylab = 'Lag Waktu (hari)', 
     zlab = 'Gamma')

Fitting cross-variogram variabel konsentrasi PM2.5 dengan kecepatan angin menghasilkan model sum-metric sebagai model dengan MSE terendah dibandingkan dengan model lainnya. Model tersebut menggunakan fungsi Gaussian pada komponen spasial, temporal, dan joint. Gambar tersebut menunjukkan perbandingan cross-variogram empiris variabel konsentrasi PM2.5 dan kecepatan angin dengan cross-variogram teoretis model sum-metric.

Predict The Test Data

sumMetric = vgmST("sumMetric", space = vgm(psill = param_summetric[1], "Gau", range = param_summetric[2], nugget = param_summetric[3]),
                  time = vgm(psill = param_summetric[4], "Mat", range = param_summetric[5], nugget = param_summetric[6]),
                  joint = vgm(psill = param_summetric[7], "Sph", range = param_summetric[8], nugget = param_summetric[9]), stAni = param_summetric[10],
                  temporalUnit = "days")
test.krige.ST = krigeST(PM25 ~ Presipitasi + Kelembaban + Angin + NO2, 
                        data = stfdf.train, 
                        newdata = stfdf.test, 
                        modelList = sumMetric)
MAPE_test = mean(abs((test.krige.ST@data$var1.pred-stfdf.test@data$PM25)/stfdf.test@data$PM25))*100
RMSE_test = sqrt(mean((test.krige.ST@data$var1.pred-stfdf.test@data$PM25)^2))

prediksi_vs_aktual = data.frame(Prediksi = test.krige.ST@data$var1.pred, Aktual = stfdf.test@data$PM25)
e = prediksi_vs_aktual$Prediksi-prediksi_vs_aktual$Aktual
#write.xlsx(prediksi_vs_aktual, "D:/Tesis/Hasil/excels/Prediksi vs Aktual Data Test.xlsx")

cat("MAPE Test: ", MAPE_test)

## MAPE Test:  0.6588764

cat("RMSE Test: ", RMSE_test)

## RMSE Test:  0.3843981

Evaluasi dilakukan menggunakan dua ukuran kekeliruan, yaitu nilai RMSE dan MAPE. Model variogram teoretis yang dipilih memberikan prediksi yang akurat pada data uji, dengan nilai MAPE sebesar 0,659% dan RMSE sebesar 0,384.

Estimasi Konsentrasi \(PM_{2.5}\) di Seluruh Wilayah Jakarta

Selanjutnya, estimasi dilakukan menggunakan model variogram dengan parameter yang telah diestimasi. Estimasi dilakukan untuk 646 grid berukuran 1 km × 1 km yang terletak secara beraturan di wilayah DKI Jakarta periode bulan Januari hingga Desember 2023.

idn.gadm = "D:/TESIS/DATA/gadm41_IDN.gpkg"
tesis.sf = st_read(idn.gadm, layer = "ADM_ADM_2")

## Reading layer `ADM_ADM_2' from data source `D:\TESIS\DATA\gadm41_IDN.gpkg' using driver `GPKG'
## Simple feature collection with 502 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.00971 ymin: -11.00761 xmax: 141.0194 ymax: 6.076941
## Geodetic CRS:  WGS 84

tesis.sf = st_transform(tesis.sf, utm48s)
wilayah = tesis.sf[tesis.sf$NAME_1 == "Jakarta Raya", ]["NAME_2"]
tesis.poly = as(wilayah[-6,], "Spatial")
tesis.poly@bbox[2, 1] <- 9290000  # Ganti `new_ymin` dengan nilai y min baru
grd = SpatialPixels(SpatialPoints(makegrid(tesis.poly, cellsize = 1000)),
                    proj4string = proj4string(tesis.poly))

n = 12 # jumlah grid waktu/periode waktu
library(xts)
tgrd = seq(min(index(stfdf.train)), max(index(stfdf.train)), length = n)

# MENGISI GRID PREDIKSI DENGAN VERIABEL PREDIKTOR MENGGUNAKAN COKRIGING SPASIAL
# UNTUK SETIAP TITIK WAKTU
bulan = c("2023-01-01", "2023-02-01", "2023-03-01", "2023-04-01", "2023-05-01", 
          "2023-06-01", "2023-07-01", "2023-08-01", "2023-09-01", "2023-10-01", 
          "2023-11-01", "2023-12-01")
prediktor_grd = c()
for(i in bulan){
  data <- stfdf.train[, i]
  gstat(id = "ch", formula = Presipitasi~1, data = data) |>
    gstat("rh", Kelembaban~1, data) |>
    gstat("no2", NO2~1, data) |>
    gstat("ws", Angin~1, data) |>
    gstat(model=vgm(psill = var(data$Presipitasi), model = "Sph", range = 30000), fill.all=T) -> g
  v <- variogram(g)
  v.fit = fit.lmc(v, g)
  plot(v, model = v.fit)
  pred <- predict(v.fit, newdata = grd)
  grd_pred = cbind(Presipitasi = pred@data$ch.pred,
                   Kelembaban = pred@data$rh.pred,
                   Angin = pred@data$ws.pred,
                   NO2 = pred@data$no2.pred)
  prediktor_grd = rbind(prediktor_grd, grd_pred)
}

## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]
## Intrinsic Correlation found. Good.
## [using ordinary cokriging]

head(prediktor_grd)

##      Presipitasi Kelembaban    Angin       NO2
## [1,]    176.0560   80.88047 4.741267 10.343881
## [2,]    176.6280   80.96468 4.734693 10.239950
## [3,]    177.2510   81.05278 4.727209 10.135030
## [4,]    177.8543   81.13615 4.719810 10.036452
## [5,]    178.4162   81.21243 4.712812  9.946155
## [6,]    178.9286   81.28097 4.706370  9.864932

# STFDF untuk new_data to predict
pred.grd = STFDF(sp = grd, time = tgrd, 
                 data = as.data.frame(prediktor_grd))

# ==============================================================PREDIKSI FULL
tesis.krige.ST = krigeST(PM25 ~ Presipitasi + Kelembaban + Angin + NO2, 
                         data = stfdf.train, newdata = pred.grd,
                         modelList = sumMetric)
colnames(tesis.krige.ST@data) = 'PM25'

min_scale = min(tesis.krige.ST@data$PM25)
max_scale = max(tesis.krige.ST@data$PM25)

# MENYIMPAN HASIL
#hasil_semua = as.data.frame(tesis.krige.ST)
#write.xlsx(hasil_semua, "D:/TESIS/Hasil/excels/Prediksi semua grid.xlsx")

# =============================================================PLOT SATU SATU
# Konversi tesis.krige.ST ke RasterLayer
raster_data <- raster(tesis.krige.ST[,2])

# Lakukan masking agar raster hanya di dalam batas wilayah
masked_raster <- mask(raster_data, tesis.poly)

# Konversi raster ke data frame
raster_df <- as.data.frame(rasterToPoints(masked_raster), xy = TRUE)
colnames(raster_df) <- c("x", "y", "PM2.5")  # Sesuaikan nama kolom

# Konversi batas wilayah (tesis.poly) ke objek sf
tesis_poly_sf <- st_as_sf(tesis.poly)
tesis_poly_sf[,"NAME_ID"] = c('Jakarta Barat', 'Jakarta Pusat', 'Jakarta Selatan',
                              'Jakarta Timur', 'Jakarta Utara')
tesis_poly_sf[,"NAME_EN"] = c('West Jakarta', 'Central Jakarta', 'South Jakarta',
                              'East Jakarta', 'North Jakarta')

# ==================================================== PLOT GABUNGAN 12 BULAN
# Vektor nama bulan dari Januari hingga Desember
#bulan <- c("January", "February", "March", "April", "May", "June",  
#           "July", "August", "September", "October", "November", "December")
bulan = c("Januari", "Februari", "Maret", "April", "Mei", "Juni", "Juli", 
          "Agustus", "September", "Oktober", "November", "Desember")
# Gabungkan semua data raster ke dalam satu data frame
all_raster_df <- do.call(rbind, lapply(1:12, function(i) {
  # Konversi tesis.krige.ST ke RasterLayer pada indeks tertentu
  raster_data <- raster(tesis.krige.ST[, i])
  
  # Lakukan masking agar raster hanya di dalam batas wilayah
  masked_raster <- mask(raster_data, tesis.poly)
  
  # Konversi raster ke data frame
  raster_df <- as.data.frame(rasterToPoints(masked_raster), xy = TRUE)
  colnames(raster_df) <- c("x", "y", "PM2.5")
  
  # Tambahkan kolom bulan
  raster_df$bulan <- bulan[i]
  raster_df$time <- i
  
  return(raster_df)
}))

# Ubah kolom 'bulan' menjadi faktor dengan urutan level yang benar
all_raster_df$bulan <- factor(all_raster_df$bulan, 
                              levels = bulan)

# Plot dengan ggplot
peta = ggplot() + 
  geom_tile(data = all_raster_df, aes(x = x, y = y, fill = PM2.5)) + 
  geom_contour(data = all_raster_df, aes(x = x, y = y, z = PM2.5, colour = stat(level)), size = 0.2, colour = "red", alpha = 0.5) +
  scale_fill_gradientn(colours = colorRampPalette(c("green4", "green", "yellow", "orangered", "red3"))(200)) +
  facet_wrap(~ bulan, ncol = 4) + 
  theme_bw() + 
  ylab("") + 
  xlab("") + 
  geom_sf(data = tesis_poly_sf, fill = NA, size=0.01, colour = "black") +  
  geom_sf_text(data = tesis_poly_sf, aes(geometry = st_centroid(geometry), label = NAME_ID),
               color = "black", size = 1.7, fontface = "bold", vjust = 1.5) + # Tambahkan label
  guides(fill = guide_colorbar(title.position = "top", title.vjust = 1,  
                               frame.colour = "black",
                               barwidth = 1,
                               barheight = 10)) +
  theme(aspect.ratio = 1, legend.position = "right", 
        text = element_text(), 
        legend.direction = "vertical",
        panel.background = element_rect(fill = "white", color = NA), 
        panel.grid.major = element_line(color = "grey95", linetype = "solid"),
        panel.grid.minor = element_line(color = "grey95", linetype = "solid"),
        axis.text.x = element_blank(), #remove x axis labels
        axis.ticks.x = element_blank(), #remove x axis ticks
        axis.text.y = element_blank(),  #remove y axis labels
        axis.ticks.y = element_blank(),  #remove y axis ticks
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7),
        strip.text = element_text(size = 7)
  ) + 
  labs(fill = "PM2.5")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

peta

## Warning: `stat(level)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Distribusi secara spasial menunjukkan bahwa konsentrasi \(PM_{2.5}\) cenderung lebih rendah pada awal tahun (Januari hingga April) dan akhir tahun (November-Desember). Konsentrasi \(PM_{2.5}\) meningkat secara signifikan pada pertengahan tahun, yaitu selama periode Mei hingga Oktober. Pola ini berkaitan dengan faktor meteorologis, seperti meningkatnya curah hujan selama musim hujan (di akhir dan awal tahun) yang melarutkan dan menghilangkan polutan, serta curah hujan yang lebih rendah selama musim kemarau, yang membatasi proses pemurnian udara. Selain itu, kecepatan angin yang lebih rendah selama musim kemarau berkontribusi pada akumulasi polutan di area tertentu.

Hasil estimasi menunjukkan bahwa konsentrasi \(PM_{2.5}\) pada tahun 2023 berada pada rentang 15,41 µg/m3 hingga 65,88 µg/m3. Indonesia menetapkan nilai ambang batas konsentrasi PM2.5 maksimal 50 µg/m3. Wilayah-wilayah dengan konsentrasi yang melebihi nilai ambang batas sehat ditunjukkan dengan warna jingga hingga merah. Terlihat bahwa konsentrasi \(PM_{2.5}\) melebihi ambang batas pada periode-periode musim kemarau.

Konsentrasi \(PM_{2.5}\) yang tinggi umumnya teramati di Jakarta Utara, Jakarta Pusat, dan sebagian wilayah Jakarta Selatan. Di Jakarta Utara dan Jakarta Pusat, peningkatan konsentrasi \(PM_{2.5}\) terjadi karena kepadatan penduduk tinggi, aktivitas industri yang tinggi, dan emisi gas dari kendaraan bermotor. Faktor-faktor tersebut juga meningkatkan konsentrasi NO2 yang berkontribusi terhadap peningkatan \(PM_{2.5}\).

Sementara itu, tingginya konsentrasi \(PM_{2.5}\) di Jakarta Selatan dipengaruhi oleh sumber lokal dan regional. Wilayah ini berbatasan langsung dengan Bogor dan Depok. Kedua wilayah tersebut memiliki tingkat urbanisasi yang tinggi dan aktivitas transportasi yang intensif. Hal ini memungkinkan adanya polutan yang bergerak ke wilayah Jakarta Selatan. Fenomena ini terutama terlihat selama musim kemarau, yaitu ketika kondisi atmosfer (kelembaban, kecepatan angin, dan curah hujan) mendukung adanya penyebaran polutan dari wilayah sekitar.